rm(list = ls())
require(stats4)

#multinomial coefficient
multinom=function(vec){
return(sum(log(1:max(1,sum(vec))))-sum(log(c(1:max(1,vec[1]),1:max(1,vec[2])))))
}



data2old=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/infoexp2_data/responses.csv", header=TRUE)
data13=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_july25_type13clean.csv", header=TRUE)
data15=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_August12_type15(clean).csv", header=TRUE)


dataframe2old=data.frame(data2old)
dataframe13=data.frame(data13)
dataframe15=data.frame(data15)

r=10

priors=c(0.5,0.6,0.75,0.85)

expframe=rbind(dataframe13,dataframe15,dataframe2old)
expframe=subset(expframe,uid>200)

#note that player 234 did not finish and so has no real data
expframe=subset(expframe,uid!=234)

#definitions
expframe$priorA=(expframe$qid==1)*0.5+(expframe$qid==5)*0.6+(expframe$qid==6)*0.75+(expframe$qid==7)*0.85
expframe$priorB=1-expframe$priorA
expframe$Aresponse=(expframe$response==6)
expframe$Bresponse=(expframe$response==7)
expframe$stateA=(expframe$state==4)
expframe$stateB=(expframe$state==10)
expframe$correct=expframe$stateA*expframe$Aresponse+expframe$stateB*expframe$Bresponse
uidlist=unique(expframe$uid)

indioutframe=data.frame(uid=vector('numeric'),likelin=vector('numeric'),errlin=vector('numeric'),likepow=vector('numeric'),errpow=vector('numeric'),likegen=vector('numeric'),errgen=vector('numeric'),likeinat=vector('numeric'),errinat=vector('numeric'))
for(indii in 1:length(uidlist)){
indiframe=subset(expframe, uid==uidlist[indii])

#Observed move frequencies

ArespAstate=vector('numeric')
BrespAstate=vector('numeric')
ArespBstate=vector('numeric')
BrespBstate=vector('numeric')

for(i in 1:4){
playframe=subset(indiframe, priorA==priors[i])
ArespAstate[i]=sum(playframe$Aresponse*playframe$stateA)
BrespAstate[i]=sum(playframe$Bresponse*playframe$stateA)
ArespBstate[i]=sum(playframe$Aresponse*playframe$stateB)
BrespBstate[i]=sum(playframe$Bresponse*playframe$stateB)
}

conmat=matrix(c(1,0, 0,1,  -1,0, 0,-1), nrow=4, byrow=TRUE)
convec=c(0.001,0.001,-0.99,-0.99 )

multinomvec=vector("numeric")
for(i in 1:4){
playframe=subset(indiframe, priorA==priors[i])
multinomvec[i]=(multinom(c(ArespAstate[i],BrespAstate[i])))+(multinom(c(ArespBstate[i],BrespBstate[i])))
}
multinomconst=sum(multinomvec)

#########################################################
#Shannon
choiceprobslin=function(kappa,treatment){
p=priors[treatment]

#A given A, A given B
cost=function(cprobs){
overallcost=kappa*(p*cprobs[1]*log(p*cprobs[1]/(p*cprobs[1]+(1-p)*cprobs[2]))+(1-p)*cprobs[2]*log((1-p)*cprobs[2]/(p*cprobs[1]+(1-p)*cprobs[2]))+p*(1-cprobs[1])*log(p*(1-cprobs[1])/(1-p*cprobs[1]-(1-p)*cprobs[2]))+(1-p)*(1-cprobs[2])*log((1-p)*(1-cprobs[2])/(1-p*cprobs[1]-(1-p)*cprobs[2])))
return(overallcost)
}

#Optimization target
negpayoff=function(cprobs){
return(-1*(r*(p*cprobs[1]+(1-p)*(1-cprobs[2])) -cost(cprobs )))
}

#Find optimal accuracies
optimout=constrOptim(rep(0.5,2),negpayoff, grad=NULL,ui=conmat,ci=convec)
cprobvec=optimout$par

return(cprobvec)

}

likelihoodlin=function(kappa,err){
cmat=matrix(c(choiceprobslin(kappa,1),choiceprobslin(kappa,2),choiceprobslin(kappa,3),choiceprobslin(kappa,4)),ncol=4)
like=multinomconst+sum(ArespAstate*log(cmat[1,]*(1-err)+err/2)+BrespAstate*log((1-cmat[1,])*(1-err)+err/2)+ArespBstate*log(cmat[2,]*(1-err)+err/2)+BrespBstate*log((1-cmat[2,])*(1-err)+err/2))
return(as.numeric(-like))
}

if(FALSE){
modellin=mle(likelihoodlin,start=list(kappa=1,err=0.001),fixed=list(err=0.001), method="L-BFGS-B",lower=c(0.001),upper=c(10000))
}else{
modellin=mle(likelihoodlin,start=list(kappa=7,err=0.1), method="L-BFGS-B",lower=c(0.001,0.001),upper=c(10000,0.999))
}

kappalin=coef(modellin)[1]
errlin=coef(modellin)[2]

lllin=logLik(modellin)[1]


##########################################
#Power Mutual Information

choiceprobspow=function(kappa,rho,treatment){
p=priors[treatment]

#A given A, A given B
cost=function(cprobs){
overallcost=(p*cprobs[1]*log(p*cprobs[1]/(p*cprobs[1]+(1-p)*cprobs[2]))+(1-p)*cprobs[2]*log((1-p)*cprobs[2]/(p*cprobs[1]+(1-p)*cprobs[2]))+p*(1-cprobs[1])*log(p*(1-cprobs[1])/(1-p*cprobs[1]-(1-p)*cprobs[2]))+(1-p)*(1-cprobs[2])*log((1-p)*(1-cprobs[2])/(1-p*cprobs[1]-(1-p)*cprobs[2])))-(p*log(p)+(1-p)*log(1-p))
return(kappa*overallcost^rho)
}

#Optimization target
negpayoff=function(cprobs){
return(-1*(r*(p*cprobs[1]+(1-p)*(1-cprobs[2])) -cost(cprobs)))
}


#Find optimal accuracies
optimout=constrOptim(c(0.6,0.5),negpayoff, grad=NULL,ui=conmat,ci=convec)
cprobvec=optimout$par

return(cprobvec)

}

likelihoodpow=function(kappa,rho,err){
cmat=matrix(c(choiceprobspow(kappa,rho,1),choiceprobspow(kappa,rho,2),choiceprobspow(kappa,rho,3),choiceprobspow(kappa,rho,4)),ncol=4)
like=multinomconst+sum(ArespAstate*log(cmat[1,]*(1-err)+err/2)+BrespAstate*log((1-cmat[1,])*(1-err)+err/2)+ArespBstate*log(cmat[2,]*(1-err)+err/2)+BrespBstate*log((1-cmat[2,])*(1-err)+err/2))
return(as.numeric(-like))
}

if(indii %in% c(8,12,20,21,42,54)){
modelpow=mle(likelihoodpow,start=list(kappa=10,rho=1,err=0.7),fixed=list(rho=4), method="L-BFGS-B",lower=c(0.001,0.01,0.001),upper=c(10000,100,0.999))
}else{
if(indii %in% c(16,31,40)){
modelpow=mle(likelihoodpow,start=list(kappa=10,rho=1,err=0.1), method="L-BFGS-B",lower=c(0.001,0.01,0.001),upper=c(10000,100,0.999))
}else{
modelpow=mle(likelihoodpow,start=list(kappa=100,rho=4,err=0.1), method="L-BFGS-B",lower=c(0.001,0.01,0.001),upper=c(10000,100,0.999))
}
}

kappapow=coef(modelpow)[1]
errpow=coef(modelpow)[3]

llpow=logLik(modelpow)[1]

############################################
#General Entropy

choiceprobsgen=function(kappa,rho,treatment){
p=priors[treatment]

#A given A, A given B
cost=function(cprobs){

if(rho==1){
overallcost=kappa*(p*cprobs[1]*log(p*cprobs[1]/(p*cprobs[1]+(1-p)*cprobs[2]))+(1-p)*cprobs[2]*log((1-p)*cprobs[2]/(p*cprobs[1]+(1-p)*cprobs[2]))+p*(1-cprobs[1])*log(p*(1-cprobs[1])/(1-p*cprobs[1]-(1-p)*cprobs[2]))+(1-p)*(1-cprobs[2])*log((1-p)*(1-cprobs[2])/(1-p*cprobs[1]-(1-p)*cprobs[2])))
}
if(rho==2){
overallcost=kappa*(-1/2*(p*log(p*cprobs[1]/(p*cprobs[1]+(1-p)*cprobs[2]))+(1-p)*log((1-p)*cprobs[2]/(p*cprobs[1]+(1-p)*cprobs[2]))+p*log(p*(1-cprobs[1])/(1-p*cprobs[1]-(1-p)*cprobs[2]))+(1-p)*log((1-p)*(1-cprobs[2])/(1-p*cprobs[1]-(1-p)*cprobs[2])))-2*log(2))
}
if(rho!=1&rho!=2){
overallcost=kappa*(1/((rho-2)*(rho-1))*2^(1-rho)*((p*cprobs[1]+(1-p)*cprobs[2])*((p*cprobs[1]/(p*cprobs[1]+(1-p)*cprobs[2]))^(2-rho)+((1-p)*cprobs[2]/(p*cprobs[1]+(1-p)*cprobs[2]))^(2-rho))+(p*(1-cprobs[1])+(1-p)*(1-cprobs[2]))*((p*(1-cprobs[1])/(1-p*cprobs[1]-(1-p)*cprobs[2]))^(2-rho)+((1-p)*(1-cprobs[2])/(1-p*cprobs[1]-(1-p)*cprobs[2]))^(2-rho) ) )-1/((rho-2)*(rho-1))-log(2))
}

return(overallcost)
}

#Optimization target
negpayoff=function(cprobs){
return(-1*(r*(p*cprobs[1]+(1-p)*(1-cprobs[2])) -cost(cprobs )))
}

#Find optimal accuracies
optimout=constrOptim(rep(0.5,2),negpayoff, grad=NULL,ui=conmat,ci=convec)
cprobvec=optimout$par

return(cprobvec)

}


likelihoodgen=function(kappa,rho,err){
cmat=matrix(c(choiceprobsgen(kappa,rho,1),choiceprobsgen(kappa,rho,2),choiceprobsgen(kappa,rho,3),choiceprobsgen(kappa,rho,4)),ncol=4)
like=multinomconst+sum(ArespAstate*log(cmat[1,]*(1-err)+err/2)+BrespAstate*log((1-cmat[1,])*(1-err)+err/2)+ArespBstate*log(cmat[2,]*(1-err)+err/2)+BrespBstate*log((1-cmat[2,])*(1-err)+err/2))
return(as.numeric(-like))
}

if(indii %in% c(20,46)){

modelgen=mle(likelihoodgen,start=list(kappa=1,rho=10,err=0.1), method="L-BFGS-B",lower=c(0.001,0.01,0.001),upper=c(10000,100,0.999))

}else{
modelgen=mle(likelihoodgen,start=list(kappa=10,rho=1,err=0.1), method="L-BFGS-B",lower=c(0.001,0.01,0.001),upper=c(10000,100,0.999))
}

kappagen=coef(modelgen)[1]
errgen=coef(modelgen)[3]

llgen=logLik(modelgen)[1]

##################################
#Inattentive Model

likelihoodinat=function(err){
cmat=matrix(c(0.5,0.5, 1,1, 1,1, 1,1),ncol=4)

like=multinomconst+sum(ArespAstate*log(cmat[1,]*(1-err)+err/2)+BrespAstate*log((1-cmat[1,])*(1-err)+err/2)+ArespBstate*log(cmat[2,]*(1-err)+err/2)+BrespBstate*log((1-cmat[2,])*(1-err)+err/2))
return(as.numeric(-like))
}

modelinat=mle(likelihoodinat,start=list(err=0.5), method="L-BFGS-B",lower=c(0.01),upper=c(0.99))
errinat=coef(modelinat)[1]
llinat=logLik(modelinat)[1]

addframe=data.frame(uid=uidlist[indii],likelin=lllin,errlin=errlin,likepow=llpow,errpow=errpow,likegen=llgen,errgen=errgen,likeinat=llinat,errinat=errinat)

indioutframe=rbind(indioutframe,addframe)

print(indii)
}


indioutframe$AIClin=indioutframe$likelin*-2+4
indioutframe$AICpow=indioutframe$likepow*-2+6
indioutframe$AICgen=indioutframe$likegen*-2+6
indioutframe$AICinat=indioutframe$likeinat*-2+2

indioutframe$minAIC=0
for(i in 1:length(indioutframe$AIClin)){
indioutframe$minAIC[i]=min(c(indioutframe$AIClin[i],indioutframe$AICpow[i],indioutframe$AICgen[i],indioutframe$AICinat[i]))
}

indioutframe$minAIClin=(indioutframe$AIClin==indioutframe$minAIC)
indioutframe$minAICpow=(indioutframe$AICpow==indioutframe$minAIC)
indioutframe$minAICgen=(indioutframe$AICgen==indioutframe$minAIC)
indioutframe$minAICinat=(indioutframe$AICinat==indioutframe$minAIC)

sum(indioutframe$minAIClin)/length(indioutframe$AIClin)
sum(indioutframe$minAICpow)/length(indioutframe$AIClin)
sum(indioutframe$minAICgen)/length(indioutframe$AIClin)
sum(indioutframe$minAICinat)/length(indioutframe$AIClin)

write.table(indioutframe, "C:/Users/neligh/Dropbox (Chapman)/Dean Stuff/Full Data Reports/experiment_3_individuals.csv", sep = ",", col.names = T, append = F)


